Disclaimer

This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.

All COOP

Coop Data

NOAA/NWS Cooperative Observer Network

For the Yampa area:

STEAMBOAT SPRINGS CO7936 HAYDEN CO3867 WALDEN CO8756 KREMMLING CO4664 CRAIG 4 SW CO1932

For the Rio Grande area:

DEL NORTE 2E CO2184 HERMIT 7 ESE CO3951 PAGOSA SPRINGS CO6258 SILVERTON CO7656 VALLECITO DAM CO8582

Downloading from the Iowa Environmental Mesonet.

COOP_stations <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/All_COOP_stations/data_raw/COOP_stations.csv", header = TRUE)

#str(COOP_stations)
COOP_stations$Date <- ymd(COOP_stations$day)

COOP_stations_clean <- COOP_stations %>% # filter for the timeframe
  addWaterYear() %>%
  mutate(daymonth = format(as.Date(Date), "%d-%m")) %>% 
  group_by(waterYear)%>% 
  mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days")))) %>%   na.omit()

COOP_stations_clean <- COOP_stations_clean %>% 
  mutate(avg_T_c = (highc+lowc)/2) %>% 
  filter(waterYear <= 2022)


write.csv(COOP_stations_clean,"C:/Users/13074/Documents/ESS580/thesis_project/All_COOP_stations/data_clean/COOP_stations_clean.csv", row.names = FALSE)

ggplot(COOP_stations_clean, aes(x = Date, y = avg_T_c, color = station_name)) +
  geom_line() + #lwd = 2) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Daily temperature (°C)') + 
  xlab('Date')

# Just checking.

Station methodology

Two stations (Del Norte and Steamboat) will be analyzed for longer spans of time (1900 - 2022), while the remaining stations will be analyzed over the same period of record that the SNOTEL stations operate at (~1980s -2022). Most SNOTEL stations analyzed do not have reliable temperature data before WY 1987 & 1984 in the northern and southern areas respectively. For this purpose we will look at data after WY 1984.

Standard Deviation

To figure out the standard deviation for each year, I want the “residual” for each daily value.

The standard deviation will be the daily residual minus the mean of the residuals by water year, summed and squared, then divided by the number of observations minus one. The square root of the resulting value of which is thus the standard deviation for the water year.

Northern Stations

HAYDEN CO3867

Detrending Data

#average water year temperature

CO3867_yearly_wy_aver <- COOP_stations_clean %>%
  filter(station == "CO3867") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO3867_daily_wy_aver <- CO3867_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO3867_daily_wy_aver <- CO3867_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO3867_daily_wy_aver$aver_day_temp))

#str(CO3867_daily_wy_aver)
# try to show all years as means. 
CO3867_daily_wy_aver2 <- CO3867_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO3867_daily_wy_aver2$date_temp <- signif(CO3867_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO3867_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO3867_standard_dev <- CO3867_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO3867_standard_dev$residual)
## [1] -7.446232e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO3867_standard_dev_all <- CO3867_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3867_standard_dev_all <- CO3867_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 11.485417
1985 11.180535
1986 9.876778
1987 9.133878
1988 11.008601
1989 11.092140
1990 10.369513
1991 10.739895
1992 10.509582
1993 10.148621
1994 10.680028
1995 9.589489
1996 9.781918
1997 10.480829
1998 9.823480
1999 9.370938
2000 9.464960
2001 11.143631
2002 11.223865
2003 10.221418
2004 9.973747
2005 9.115467
2006 10.817760
2007 11.295537
2008 10.724288
2009 9.396883
2010 10.757495
2011 10.786836
2012 9.858540
2013 11.469037
2014 10.317811
2015 9.518729
2016 10.994901
2017 9.957091
2018 9.664345
2019 10.669019
2020 10.820406
2021 10.982860
2022 10.516975 
#CALLING THIS something different
CO3867_all_V2 <- ggplot(CO3867_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')


print(CO3867_all_V2)

CO3867_sd_mk <- mk.test(CO3867_standard_dev_all$sd_2)
print(CO3867_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all$sd_2
## z = -0.096775, n = 39, p-value = 0.9229
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##   -9.00000000 6833.66666667   -0.01214575
CO3867_sd_sens <- sens.slope(CO3867_standard_dev_all$sd_2)
print(CO3867_sd_sens)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all$sd_2
## z = -0.096775, n = 39, p-value = 0.9229
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02209965  0.01989011
## sample estimates:
##   Sen's slope 
## -0.0007800241

Summer temperature standard deviation

CO3867_standard_dev_all_summer <- CO3867_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3867_standard_dev_all_summer <- CO3867_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.618978
1985 2.833661
1986 1.841107
1987 2.301978
1988 2.094272
1989 3.019808
1990 2.776394
1991 2.268129
1992 2.990114
1993 2.697978
1994 2.543786
1995 3.295977
1996 2.378301
1997 2.201010
1998 3.847913
1999 2.605061
2000 2.576684
2001 3.026525
2002 2.742324
2003 3.677107
2004 2.830979
2005 3.705172
2006 2.325716
2007 3.226804
2008 3.564297
2009 2.757671
2010 2.749515
2011 2.736488
2012 2.756577
2013 2.504478
2014 2.697583
2015 2.091126
2016 2.396017
2017 2.530960
2018 2.829001
2019 3.318852
2020 2.939140
2021 3.007404
2022 2.815852 
ggplot(CO3867_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO3867_sd_mk_summer <- mk.test(CO3867_standard_dev_all_summer$sd_2)
print(CO3867_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_summer$sd_2
## z = 0.65323, n = 39, p-value = 0.5136
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.500000e+01 6.833667e+03 7.422402e-02
CO3867_sd_sens_summer <- sens.slope(CO3867_standard_dev_all_summer$sd_2)
print(CO3867_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_summer$sd_2
## z = 0.65323, n = 39, p-value = 0.5136
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01003205  0.01812277
## sample estimates:
## Sen's slope 
##  0.00511551

Winter temperature standard deviation

CO3867_standard_dev_all_winter <- CO3867_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3867_standard_dev_all_winter <- CO3867_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3867_standard_dev_all_winter <- CO3867_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.752371
1985 6.992625
1986 7.021863
1987 4.746588
1988 6.458530
1989 7.858157
1990 6.065225
1991 6.992195
1992 6.622718
1993 5.995240
1994 5.966573
1995 5.927126
1996 5.644068
1997 6.264780
1998 4.825588
1999 6.279732
2000 5.566801
2001 5.236716
2002 6.735051
2003 5.328239
2004 6.396353
2005 4.989144
2006 6.367099
2007 7.523509
2008 6.357499
2009 5.844773
2010 6.267521
2011 7.004881
2012 5.213377
2013 7.479074
2014 6.512199
2015 6.697790
2016 6.577107
2017 7.276273
2018 5.600155
2019 5.454007
2020 5.743898
2021 5.775806
2022 6.062583 
ggplot(CO3867_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO3867_sd_mk_winter <- mk.test(CO3867_standard_dev_all_winter$sd_2)
print(CO3867_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_winter$sd_2
## z = -0.84678, n = 39, p-value = 0.3971
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -71.00000000 6833.66666667   -0.09581646
CO3867_sd_sens_winter <- sens.slope(CO3867_standard_dev_all_winter$sd_2)
print(CO3867_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_winter$sd_2
## z = -0.84678, n = 39, p-value = 0.3971
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.03164878  0.01862392
## sample estimates:
##  Sen's slope 
## -0.009441259

Spring and Fall temperature standard deviation

CO3867_standard_dev_all_spring <- CO3867_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3867_standard_dev_all_spring <- CO3867_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.383271
1985 4.121032
1986 3.945268
1987 4.252628
1988 4.686791
1989 4.184255
1990 3.594536
1991 4.930077
1992 3.554712
1993 5.092515
1994 4.882068
1995 3.347829
1996 4.640987
1997 6.021795
1998 4.152854
1999 5.103695
2000 4.548147
2001 4.268838
2002 3.799616
2003 5.469983
2004 3.661611
2005 4.356924
2006 4.109801
2007 4.274371
2008 5.406461
2009 5.078561
2010 4.761126
2011 3.716299
2012 3.837497
2013 6.187307
2014 5.221806
2015 3.236390
2016 3.966077
2017 4.259327
2018 4.506738
2019 3.634707
2020 5.392938
2021 4.684954
2022 4.794863 
ggplot(CO3867_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO3867_sd_mk_spring <- mk.test(CO3867_standard_dev_all_spring$sd_2)
print(CO3867_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_spring$sd_2
## z = 0.21774, n = 39, p-value = 0.8276
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.900000e+01 6.833667e+03 2.564103e-02
CO3867_sd_sens_spring <- sens.slope(CO3867_standard_dev_all_spring$sd_2)
print(CO3867_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_spring$sd_2
## z = 0.21774, n = 39, p-value = 0.8276
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02256156  0.02462188
## sample estimates:
## Sen's slope 
## 0.002681139

Fall temperature standard deviation

CO3867_standard_dev_all_fall <- CO3867_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3867_standard_dev_all_fall <- CO3867_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3867_standard_dev_all_fall <- CO3867_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3867_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.261658
1985 6.543494
1986 4.185360
1987 4.088324
1988 4.045865
1989 3.385108
1990 5.695730
1991 4.597785
1992 5.560309
1993 3.458722
1994 5.412967
1995 5.437667
1996 5.191045
1997 7.215270
1998 6.225676
1999 4.095995
2000 4.936885
2001 5.012137
2002 4.507905
2003 4.960323
2004 4.671952
2005 4.452681
2006 4.089773
2007 5.747986
2008 4.077807
2009 4.885493
2010 5.812062
2011 4.255662
2012 4.967710
2013 5.812048
2014 5.746594
2015 3.902345
2016 3.915061
2017 4.601843
2018 5.532722
2019 5.699879
2020 8.047056
2021 7.309303
2022 6.024886 
ggplot(CO3867_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3867 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO3867_sd_mk_fall <- mk.test(CO3867_standard_dev_all_fall$sd_2)
print(CO3867_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO3867_standard_dev_all_fall$sd_2
## z = 1.5242, n = 39, p-value = 0.1275
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##          S       varS        tau 
##  127.00000 6833.66667    0.17139
CO3867_sd_sens_fall <- sens.slope(CO3867_standard_dev_all_fall$sd_2)
print(CO3867_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO3867_standard_dev_all_fall$sd_2
## z = 1.5242, n = 39, p-value = 0.1275
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.006642115  0.055193772
## sample estimates:
## Sen's slope 
##  0.02212322

WALDEN CO8756

Detrending Data

#average water year temperature

CO8756_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO8756") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO8756_daily_wy_aver <- CO8756_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO8756_daily_wy_aver <- CO8756_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO8756_daily_wy_aver$aver_day_temp))

#str(CO8756_daily_wy_aver)
# try to show all years as means. 
CO8756_daily_wy_aver2 <- CO8756_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO8756_daily_wy_aver2$date_temp <- signif(CO8756_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO8756_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO8756_standard_dev <- CO8756_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO8756_standard_dev$residual)
## [1] 1.497597e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO8756_standard_dev_all <- CO8756_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8756_standard_dev_all <- CO8756_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.274351
1985 9.525998
1986 8.779071
1987 8.487845
1988 10.339904
1989 9.651614
1990 9.041034
1991 9.094379
1992 8.889496
1993 9.021548
1994 9.472534
1995 8.585995
1996 9.196590
1997 9.540266
1998 9.195245
1999 8.607493
2000 8.580891
2001 9.998693
2002 9.568990
2003 9.169788
2004 9.135617
2005 9.362395
2006 9.918959
2007 10.583525
2008 10.168525
2009 8.515790
2010 10.879203
2011 10.106370
2012 9.301494
2013 10.654885
2014 9.305167
2015 8.952216
2016 10.068499
2017 9.421729
2018 8.417010
2019 9.947851
2020 10.248075
2021 9.934583
2022 10.543226 
CO8756_all <- ggplot(CO8756_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756_all

CO8756_sd_mk <- mk.test(CO8756_standard_dev_all$sd_2)
print(CO8756_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all$sd_2
## z = 1.8629, n = 39, p-value = 0.06247
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  155.0000000 6833.6666667    0.2091768
CO8756_sd_sens <- sens.slope(CO8756_standard_dev_all$sd_2)
print(CO8756_sd_sens)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all$sd_2
## z = 1.8629, n = 39, p-value = 0.06247
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.002285023  0.041751906
## sample estimates:
## Sen's slope 
##  0.01980816

Summer temperature standard deviation

CO8756_standard_dev_all_summer <- CO8756_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8756_standard_dev_all_summer <- CO8756_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.098045
1985 2.400143
1986 1.746046
1987 2.476687
1988 2.215917
1989 2.831668
1990 2.497958
1991 1.830750
1992 2.688531
1993 2.466431
1994 2.130837
1995 2.887343
1996 2.187068
1997 1.919971
1998 3.590781
1999 2.562795
2000 2.498102
2001 3.070289
2002 2.699482
2003 3.257005
2004 2.445442
2005 3.331090
2006 2.618565
2007 3.150664
2008 3.275320
2009 2.537241
2010 2.571604
2011 2.584250
2012 2.702433
2013 2.595520
2014 2.526453
2015 2.250917
2016 2.303665
2017 2.443964
2018 2.668405
2019 3.050127
2020 2.352849
2021 2.765161
2022 3.145041 
ggplot(CO8756_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO8756_sd_mk_summer <- mk.test(CO8756_standard_dev_all_summer$sd_2)
print(CO8756_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_summer$sd_2
## z = 1.4032, n = 39, p-value = 0.1605
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  117.0000000 6833.6666667    0.1578947
CO8756_sd_sens_summer <- sens.slope(CO8756_standard_dev_all_summer$sd_2)
print(CO8756_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_summer$sd_2
## z = 1.4032, n = 39, p-value = 0.1605
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.003752647  0.018993274
## sample estimates:
## Sen's slope 
## 0.007350181

Winter temperature standard deviation

CO8756_standard_dev_all_winter <- CO8756_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8756_standard_dev_all_winter <- CO8756_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8756_standard_dev_all_winter <- CO8756_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.388216
1985 6.145140
1986 6.192094
1987 4.775302
1988 6.147741
1989 7.299660
1990 5.674100
1991 5.806674
1992 5.383528
1993 5.497675
1994 5.335567
1995 5.625263
1996 6.227290
1997 6.355164
1998 4.852754
1999 5.892465
2000 5.410936
2001 5.139805
2002 5.879976
2003 5.478155
2004 6.577846
2005 5.799475
2006 6.730566
2007 8.012348
2008 6.720676
2009 5.728422
2010 7.020769
2011 7.094309
2012 5.681759
2013 7.926274
2014 5.973882
2015 7.409992
2016 7.063565
2017 7.624513
2018 5.533485
2019 5.538138
2020 6.860450
2021 5.784767
2022 7.619323 
ggplot(CO8756_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO8756_sd_mk_winter <- mk.test(CO8756_standard_dev_all_winter$sd_2)
print(CO8756_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_winter$sd_2
## z = 2.129, n = 39, p-value = 0.03325
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  177.0000000 6833.6666667    0.2388664
CO8756_sd_sens_winter <- sens.slope(CO8756_standard_dev_all_winter$sd_2)
print(CO8756_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_winter$sd_2
## z = 2.129, n = 39, p-value = 0.03325
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.003425964 0.056677063
## sample estimates:
## Sen's slope 
##  0.02683425

Spring and Fall temperature standard deviation

CO8756_standard_dev_all_spring <- CO8756_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8756_standard_dev_all_spring <- CO8756_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.277972
1985 3.643710
1986 3.577343
1987 4.415793
1988 4.738199
1989 4.303637
1990 3.485296
1991 5.087750
1992 3.577874
1993 4.561083
1994 4.599138
1995 3.616581
1996 4.476834
1997 5.644753
1998 3.882234
1999 5.170280
2000 4.607206
2001 4.220929
2002 3.486062
2003 5.064096
2004 3.950975
2005 4.795992
2006 4.043212
2007 4.573008
2008 6.130959
2009 5.842557
2010 4.970651
2011 4.278117
2012 3.898188
2013 6.407499
2014 5.727336
2015 3.739759
2016 4.626576
2017 4.679278
2018 4.219974
2019 3.872194
2020 6.157083
2021 4.799917
2022 5.658442 
ggplot(CO8756_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO8756_sd_mk_spring <- mk.test(CO8756_standard_dev_all_spring$sd_2)
print(CO8756_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_spring$sd_2
## z = 1.6694, n = 39, p-value = 0.09504
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  139.0000000 6833.6666667    0.1875843
CO8756_sd_sens_spring <- sens.slope(CO8756_standard_dev_all_spring$sd_2)
print(CO8756_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_spring$sd_2
## z = 1.6694, n = 39, p-value = 0.09504
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.005252066  0.044056161
## sample estimates:
## Sen's slope 
##   0.0192142

Fall temperature standard deviation

CO8756_standard_dev_all_fall <- CO8756_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8756_standard_dev_all_fall <- CO8756_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8756_standard_dev_all_fall <- CO8756_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8756_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.350339
1985 6.081834
1986 3.710357
1987 3.676327
1988 3.618004
1989 2.969292
1990 5.595828
1991 4.947550
1992 5.118606
1993 3.362972
1994 5.103288
1995 5.083007
1996 4.841146
1997 6.376457
1998 6.010537
1999 3.800887
2000 4.775100
2001 4.400595
2002 4.326695
2003 4.410416
2004 3.882484
2005 3.835660
2006 4.037489
2007 6.196685
2008 4.173587
2009 5.353327
2010 5.720666
2011 4.519050
2012 5.374169
2013 6.193739
2014 5.535980
2015 3.910185
2016 3.996059
2017 4.415269
2018 5.885926
2019 5.792724
2020 7.972408
2021 7.850006
2022 5.133755 
ggplot(CO8756_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8756 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO8756_sd_mk_fall <- mk.test(CO8756_standard_dev_all_fall$sd_2)
print(CO8756_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO8756_standard_dev_all_fall$sd_2
## z = 2.2258, n = 39, p-value = 0.02603
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  185.0000000 6833.6666667    0.2496626
CO8756_sd_sens_fall <- sens.slope(CO8756_standard_dev_all_fall$sd_2)
print(CO8756_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO8756_standard_dev_all_fall$sd_2
## z = 2.2258, n = 39, p-value = 0.02603
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.004910026 0.068048874
## sample estimates:
## Sen's slope 
##  0.03490952

KREMMLING CO4664

Detrending Data

#average water year temperature

CO4664_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO4664") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO4664_daily_wy_aver <- CO4664_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO4664_daily_wy_aver <- CO4664_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO4664_daily_wy_aver$aver_day_temp))

#str(CO4664_daily_wy_aver)
# try to show all years as means. 
CO4664_daily_wy_aver2 <- CO4664_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO4664_daily_wy_aver2$date_temp <- signif(CO4664_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO4664_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO4664_standard_dev <- CO4664_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO4664_standard_dev$residual)
## [1] -1.49728e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO4664_standard_dev_all <- CO4664_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO4664_standard_dev_all <- CO4664_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 11.842808
1985 10.534290
1986 10.076518
1987 9.278036
1988 11.625055
1989 9.347699
1990 8.896884
1991 11.059919
1992 10.593674
1993 11.082078
1994 10.858112
1995 10.294300
1996 10.772898
1997 10.417175
1998 10.004330
1999 9.816385
2000 9.827327
2001 11.308940
2002 11.334081
2003 10.241340
2004 9.836909
2005 10.949455
2006 11.325117
2007 11.322739
2008 11.113829
2009 9.570323
2010 11.562640
2011 10.555622
2012 10.157042
2013 12.623254
2014 10.708303
2015 9.799999
2016 11.086904
2017 9.803356
2018 9.386863
2019 10.644464
2020 11.277542
2021 10.570883
2022 12.052823 
ggplot(CO4664_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664_sd_mk <- mk.test(CO4664_standard_dev_all$sd_2)
print(CO4664_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all$sd_2
## z = 0.65323, n = 39, p-value = 0.5136
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.500000e+01 6.833667e+03 7.422402e-02
CO4664_sd_sens <- sens.slope(CO4664_standard_dev_all$sd_2)
print(CO4664_sd_sens)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all$sd_2
## z = 0.65323, n = 39, p-value = 0.5136
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01825695  0.03598021
## sample estimates:
## Sen's slope 
##  0.01033134

Summer temperature standard deviation

CO4664_standard_dev_all_summer <- CO4664_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO4664_standard_dev_all_summer <- CO4664_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.215574
1985 2.475205
1986 2.101595
1987 2.347529
1988 2.200798
1989 2.800050
1990 2.226887
1991 2.076735
1992 2.867749
1993 2.776034
1994 2.281763
1995 3.141926
1996 2.314091
1997 2.451858
1998 3.543159
1999 2.878018
2000 2.506292
2001 3.257659
2002 2.660263
2003 3.268403
2004 3.015509
2005 3.361711
2006 2.706764
2007 3.234749
2008 3.284052
2009 2.592534
2010 2.867244
2011 2.701616
2012 2.305239
2013 2.744107
2014 2.749757
2015 2.245707
2016 2.532862
2017 2.363437
2018 2.522825
2019 3.257400
2020 2.507018
2021 2.477895
2022 2.912039 
ggplot(CO4664_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO4664_sd_mk_summer <- mk.test(CO4664_standard_dev_all_summer$sd_2)
print(CO4664_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_summer$sd_2
## z = 0.87097, n = 39, p-value = 0.3838
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 7.300000e+01 6.833667e+03 9.851552e-02
CO4664_sd_sens_summer <- sens.slope(CO4664_standard_dev_all_summer$sd_2)
print(CO4664_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_summer$sd_2
## z = 0.87097, n = 39, p-value = 0.3838
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.008550703  0.017390398
## sample estimates:
## Sen's slope 
## 0.004689828

Winter temperature standard deviation

CO4664_standard_dev_all_winter <- CO4664_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO4664_standard_dev_all_winter <- CO4664_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO4664_standard_dev_all_winter <- CO4664_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 8.415912
1985 6.582282
1986 7.474738
1987 5.125150
1988 7.328129
1989 6.353786
1990 4.920860
1991 7.488192
1992 6.831213
1993 6.969575
1994 6.705696
1995 7.240966
1996 7.611583
1997 7.031740
1998 5.434719
1999 7.362196
2000 5.892683
2001 6.231018
2002 7.028370
2003 5.938219
2004 7.200707
2005 7.097081
2006 7.531862
2007 7.995449
2008 6.745011
2009 6.431922
2010 7.100432
2011 7.036948
2012 5.800580
2013 9.665351
2014 7.332683
2015 6.994797
2016 7.667858
2017 7.713234
2018 5.787009
2019 5.846507
2020 7.142392
2021 5.555605
2022 8.574035 
ggplot(CO4664_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO4664_sd_mk_winter <- mk.test(CO4664_standard_dev_all_winter$sd_2)
print(CO4664_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_winter$sd_2
## z = 0.55646, n = 39, p-value = 0.5779
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   47.0000000 6833.6666667    0.0634278
CO4664_sd_sens_winter <- sens.slope(CO4664_standard_dev_all_winter$sd_2)
print(CO4664_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_winter$sd_2
## z = 0.55646, n = 39, p-value = 0.5779
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02129328  0.03407467
## sample estimates:
## Sen's slope 
## 0.006437321

Spring and Fall temperature standard deviation

CO4664_standard_dev_all_spring <- CO4664_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO4664_standard_dev_all_spring <- CO4664_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.664550
1985 3.571570
1986 3.633420
1987 3.899034
1988 4.324779
1989 4.581225
1990 3.296714
1991 4.979723
1992 3.459486
1993 4.575851
1994 4.551999
1995 3.446808
1996 4.689661
1997 5.177740
1998 4.305844
1999 5.382816
2000 4.844386
2001 4.507415
2002 3.565434
2003 5.016175
2004 3.880826
2005 4.419094
2006 4.085707
2007 4.120342
2008 5.147569
2009 5.606380
2010 4.700846
2011 4.469301
2012 3.836638
2013 5.660520
2014 5.479750
2015 3.149099
2016 3.985349
2017 4.339221
2018 4.200372
2019 3.509799
2020 5.129284
2021 4.568910
2022 4.989100 
ggplot(CO4664_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO4664_sd_mk_spring <- mk.test(CO4664_standard_dev_all_spring$sd_2)
print(CO4664_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_spring$sd_2
## z = 0.7742, n = 39, p-value = 0.4388
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   65.0000000 6833.6666667    0.0877193
CO4664_sd_sens_spring <- sens.slope(CO4664_standard_dev_all_spring$sd_2)
print(CO4664_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_spring$sd_2
## z = 0.7742, n = 39, p-value = 0.4388
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01311787  0.03518940
## sample estimates:
## Sen's slope 
##  0.01059481

Fall temperature standard deviation

CO4664_standard_dev_all_fall <- CO4664_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO4664_standard_dev_all_fall <- CO4664_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO4664_standard_dev_all_fall <- CO4664_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO4664_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.800534
1985 5.716268
1986 3.740969
1987 3.729478
1988 3.543089
1989 2.956713
1990 5.580755
1991 4.544218
1992 5.419107
1993 3.661874
1994 5.209689
1995 5.283122
1996 5.399883
1997 6.298959
1998 6.262325
1999 4.344379
2000 5.423498
2001 4.444621
2002 4.685162
2003 4.620854
2004 4.021034
2005 3.868084
2006 3.716716
2007 6.109469
2008 4.351967
2009 5.241621
2010 5.836942
2011 4.464404
2012 5.751180
2013 5.656729
2014 5.287331
2015 3.708606
2016 3.807991
2017 4.343185
2018 5.508533
2019 5.380021
2020 7.334260
2021 7.658668
2022 5.204699 
ggplot(CO4664_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO4664 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO4664_sd_mk_fall <- mk.test(CO4664_standard_dev_all_fall$sd_2)
print(CO4664_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO4664_standard_dev_all_fall$sd_2
## z = 1.621, n = 39, p-value = 0.105
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  135.0000000 6833.6666667    0.1821862
CO4664_sd_sens_fall <- sens.slope(CO4664_standard_dev_all_fall$sd_2)
print(CO4664_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO4664_standard_dev_all_fall$sd_2
## z = 1.621, n = 39, p-value = 0.105
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.00290456  0.06127126
## sample estimates:
## Sen's slope 
##  0.02891896

CRAIG 4 SW CO1932

Detrending Data

#average water year temperature

CO1932_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO1932") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO1932_daily_wy_aver <- CO1932_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO1932_daily_wy_aver <- CO1932_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO1932_daily_wy_aver$aver_day_temp))

#str(CO1932_daily_wy_aver)
# try to show all years as means. 
CO1932_daily_wy_aver2 <- CO1932_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO1932_daily_wy_aver2$date_temp <- signif(CO1932_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO1932_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO1932_standard_dev <- CO1932_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO1932_standard_dev$residual)
## [1] -8.879044e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO1932_standard_dev_all <- CO1932_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO1932_standard_dev_all <- CO1932_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 11.488319
1985 11.455082
1986 9.864992
1987 9.585951
1988 11.246784
1989 11.555146
1990 10.344382
1991 10.858624
1992 10.472265
1993 10.301132
1994 10.536316
1995 9.310312
1996 9.544994
1997 10.099325
1998 9.847825
1999 9.088993
2000 9.304220
2001 11.005474
2002 11.289134
2003 9.843375
2004 10.103283
2005 8.954807
2006 10.204469
2007 11.470209
2008 10.986494
2009 9.784410
2010 10.911922
2011 10.981962
2012 10.280547
2013 11.815300
2014 10.513937
2015 9.749572
2016 11.080780
2017 10.244804
2018 9.882973
2019 10.893427
2020 11.218502
2021 11.125873
2022 11.083751 
ggplot(CO1932_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932_sd_mk <- mk.test(CO1932_standard_dev_all$sd_2)
print(CO1932_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all$sd_2
## z = 0.45968, n = 39, p-value = 0.6457
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 3.900000e+01 6.833667e+03 5.263158e-02
CO1932_sd_sens <- sens.slope(CO1932_standard_dev_all$sd_2)
print(CO1932_sd_sens)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all$sd_2
## z = 0.45968, n = 39, p-value = 0.6457
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01573645  0.03213936
## sample estimates:
## Sen's slope 
## 0.007273959

Summer temperature standard deviation

CO1932_standard_dev_all_summer <- CO1932_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO1932_standard_dev_all_summer <- CO1932_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.771767
1985 3.061372
1986 2.236917
1987 2.755026
1988 2.923342
1989 3.447128
1990 3.285824
1991 2.610959
1992 3.341025
1993 3.172880
1994 2.326110
1995 3.346459
1996 2.373395
1997 2.402126
1998 4.132837
1999 2.823431
2000 2.788161
2001 3.157886
2002 2.983087
2003 3.610345
2004 3.041601
2005 3.984741
2006 2.667119
2007 3.620785
2008 4.106893
2009 3.096088
2010 3.127090
2011 3.104357
2012 2.684607
2013 2.760181
2014 3.435432
2015 2.536963
2016 2.475548
2017 2.685470
2018 3.126233
2019 3.602189
2020 3.309112
2021 2.950159
2022 3.367256 
ggplot(CO1932_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO1932_sd_mk_summer <- mk.test(CO1932_standard_dev_all_summer$sd_2)
print(CO1932_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_summer$sd_2
## z = 0.48387, n = 39, p-value = 0.6285
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.100000e+01 6.833667e+03 5.533063e-02
CO1932_sd_sens_summer <- sens.slope(CO1932_standard_dev_all_summer$sd_2)
print(CO1932_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_summer$sd_2
## z = 0.48387, n = 39, p-value = 0.6285
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01193446  0.01886935
## sample estimates:
## Sen's slope 
## 0.003743971

Winter temperature standard deviation

CO1932_standard_dev_all_winter <- CO1932_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO1932_standard_dev_all_winter <- CO1932_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO1932_standard_dev_all_winter <- CO1932_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.763802
1985 6.756828
1986 6.641206
1987 5.189122
1988 6.710848
1989 8.371758
1990 6.190275
1991 7.344429
1992 6.468446
1993 6.185449
1994 5.749516
1995 5.899190
1996 5.778460
1997 6.082748
1998 4.827772
1999 5.919596
2000 5.395365
2001 5.047342
2002 6.813720
2003 4.981841
2004 6.365125
2005 4.949543
2006 5.965877
2007 7.722619
2008 6.497950
2009 6.262626
2010 6.318629
2011 7.345943
2012 5.404992
2013 7.814832
2014 6.603225
2015 7.017268
2016 6.895786
2017 7.501350
2018 5.655495
2019 5.395841
2020 5.893671
2021 5.718352
2022 6.929426 
ggplot(CO1932_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO1932_sd_mk_winter <- mk.test(CO1932_standard_dev_all_winter$sd_2)
print(CO1932_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_winter$sd_2
## z = -0.096775, n = 39, p-value = 0.9229
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##   -9.00000000 6833.66666667   -0.01214575
CO1932_sd_sens_winter <- sens.slope(CO1932_standard_dev_all_winter$sd_2)
print(CO1932_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_winter$sd_2
## z = -0.096775, n = 39, p-value = 0.9229
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02907046  0.02725867
## sample estimates:
##  Sen's slope 
## -0.001356466

Spring and Fall temperature standard deviation

CO1932_standard_dev_all_spring <- CO1932_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO1932_standard_dev_all_spring <- CO1932_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.582355
1985 4.318045
1986 4.388570
1987 4.465872
1988 5.099837
1989 4.890259
1990 3.471821
1991 4.812329
1992 3.789584
1993 5.391180
1994 5.063140
1995 3.114639
1996 4.628536
1997 6.362497
1998 4.065451
1999 5.287371
2000 4.653658
2001 4.693349
2002 3.742442
2003 5.293587
2004 3.745658
2005 4.598819
2006 4.480922
2007 4.442672
2008 5.838826
2009 5.893140
2010 4.949366
2011 4.358525
2012 4.213837
2013 6.736760
2014 5.609859
2015 3.328040
2016 4.175292
2017 4.244721
2018 4.759261
2019 4.013562
2020 5.826422
2021 5.376569
2022 5.094627 
ggplot(CO1932_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO1932_sd_mk_spring <- mk.test(CO1932_standard_dev_all_spring$sd_2)
print(CO1932_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_spring$sd_2
## z = 0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.700000e+01 6.833667e+03 3.643725e-02
CO1932_sd_sens_spring <- sens.slope(CO1932_standard_dev_all_spring$sd_2)
print(CO1932_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_spring$sd_2
## z = 0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02447941  0.03019881
## sample estimates:
## Sen's slope 
## 0.005866851

Fall temperature standard deviation

CO1932_standard_dev_all_fall <- CO1932_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO1932_standard_dev_all_fall <- CO1932_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO1932_standard_dev_all_fall <- CO1932_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO1932_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 4.796075
1985 7.191544
1986 4.658906
1987 4.781504
1988 4.662647
1989 3.604551
1990 6.373658
1991 5.180739
1992 6.002025
1993 3.777690
1994 5.761621
1995 5.333754
1996 5.364837
1997 7.193156
1998 6.298731
1999 3.939938
2000 5.215216
2001 4.717350
2002 4.938670
2003 4.472349
2004 4.949713
2005 3.949634
2006 4.273027
2007 6.338742
2008 4.419023
2009 5.501496
2010 6.214867
2011 5.054181
2012 5.535662
2013 6.427863
2014 5.887214
2015 4.305998
2016 4.201860
2017 4.967837
2018 6.405513
2019 6.320429
2020 8.744406
2021 8.067855
2022 5.900400 
ggplot(CO1932_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO1932 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO1932_sd_mk_fall <- mk.test(CO1932_standard_dev_all_fall$sd_2)
print(CO1932_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO1932_standard_dev_all_fall$sd_2
## z = 1.6452, n = 39, p-value = 0.09993
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  137.0000000 6833.6666667    0.1848853
CO1932_sd_sens_fall <- sens.slope(CO1932_standard_dev_all_fall$sd_2)
print(CO1932_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO1932_standard_dev_all_fall$sd_2
## z = 1.6452, n = 39, p-value = 0.09993
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.008188562  0.060716454
## sample estimates:
## Sen's slope 
##  0.02641383

STEAMBOAT SPRINGS CO7936

Detrending Data

#average water year temperature

CO7936_yearly_wy_aver <- COOP_stations_clean %>%
  filter(station == "CO7936") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO7936_daily_wy_aver <- CO7936_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO7936_daily_wy_aver <- CO7936_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO7936_daily_wy_aver$aver_day_temp))

#str(CO7936_daily_wy_aver)
# try to show all years as means. 
CO7936_daily_wy_aver2 <- CO7936_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO7936_daily_wy_aver2$date_temp <- signif(CO7936_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO7936_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO7936_standard_dev <- CO7936_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO7936_standard_dev$residual)
## [1] 1.536499e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO7936_standard_dev_all <- CO7936_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7936_standard_dev_all <- CO7936_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.500315
1985 10.331049
1986 9.546770
1987 9.102168
1988 10.958152
1989 10.863926
1990 10.023350
1991 10.418768
1992 9.629944
1993 10.268430
1994 10.466511
1995 9.467851
1996 9.705887
1997 10.248243
1998 10.461905
1999 9.873779
2000 9.769677
2001 11.360586
2002 11.435819
2003 10.114558
2004 9.833531
2005 9.123187
2006 10.735583
2007 10.818504
2008 10.266246
2009 9.317293
2010 10.270779
2011 10.245595
2012 9.836316
2013 11.024487
2014 9.892355
2015 8.911855
2016 10.397083
2017 9.675927
2018 9.395575
2019 10.325819
2020 10.362244
2021 10.332149
2022 10.356318 
#CALLING THIS something different
CO7936_all_V2 <- ggplot(CO7936_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')


print(CO7936_all_V2)

CO7936_sd_mk <- mk.test(CO7936_standard_dev_all$sd_2)
print(CO7936_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all$sd_2
## z = -0.24194, n = 39, p-value = 0.8088
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -21.00000000 6833.66666667   -0.02834008
CO7936_sd_sens <- sens.slope(CO7936_standard_dev_all$sd_2)
print(CO7936_sd_sens)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all$sd_2
## z = -0.24194, n = 39, p-value = 0.8088
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02030299  0.01505127
## sample estimates:
##  Sen's slope 
## -0.002523423

Summer temperature standard deviation

CO7936_standard_dev_all_summer <- CO7936_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7936_standard_dev_all_summer <- CO7936_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.172889
1985 2.475447
1986 1.756895
1987 2.256401
1988 2.285152
1989 3.360473
1990 2.852068
1991 2.058611
1992 3.096955
1993 3.067035
1994 2.426131
1995 3.305683
1996 2.377386
1997 2.398553
1998 3.580276
1999 2.797057
2000 2.904944
2001 3.460685
2002 3.017018
2003 3.467162
2004 3.122014
2005 3.101184
2006 2.707708
2007 3.544183
2008 3.549712
2009 2.728170
2010 2.607813
2011 2.753046
2012 2.597680
2013 2.552035
2014 3.020503
2015 2.073783
2016 2.293315
2017 2.437804
2018 2.919615
2019 3.330728
2020 2.941152
2021 2.660758
2022 3.097638 
ggplot(CO7936_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO7936_sd_mk_summer <- mk.test(CO7936_standard_dev_all_summer$sd_2)
print(CO7936_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_summer$sd_2
## z = 0.58065, n = 39, p-value = 0.5615
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.900000e+01 6.833667e+03 6.612686e-02
CO7936_sd_sens_summer <- sens.slope(CO7936_standard_dev_all_summer$sd_2)
print(CO7936_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_summer$sd_2
## z = 0.58065, n = 39, p-value = 0.5615
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01176159  0.01936204
## sample estimates:
## Sen's slope 
## 0.005147536

Winter temperature standard deviation

CO7936_standard_dev_all_winter <- CO7936_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7936_standard_dev_all_winter <- CO7936_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7936_standard_dev_all_winter <- CO7936_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.646515
1985 6.458234
1986 6.568330
1987 4.514914
1988 6.587116
1989 7.659194
1990 6.368994
1991 7.092971
1992 6.422494
1993 6.297419
1994 5.731104
1995 5.934439
1996 5.843102
1997 5.830412
1998 5.402242
1999 6.520192
2000 5.513113
2001 5.891148
2002 7.045910
2003 5.302490
2004 6.473448
2005 5.182297
2006 6.100004
2007 7.024869
2008 5.906728
2009 5.640356
2010 5.983123
2011 6.774907
2012 5.146322
2013 7.155593
2014 5.857218
2015 5.946991
2016 5.826514
2017 7.325019
2018 5.197205
2019 5.312714
2020 5.825474
2021 5.355598
2022 6.424842 
ggplot(CO7936_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO7936_sd_mk_winter <- mk.test(CO7936_standard_dev_all_winter$sd_2)
print(CO7936_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_winter$sd_2
## z = -1.7661, n = 39, p-value = 0.07737
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -147.0000000 6833.6666667   -0.1983806
CO7936_sd_sens_winter <- sens.slope(CO7936_standard_dev_all_winter$sd_2)
print(CO7936_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_winter$sd_2
## z = -1.7661, n = 39, p-value = 0.07737
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.036423474  0.002225703
## sample estimates:
## Sen's slope 
## -0.01995043

Spring and Fall temperature standard deviation

CO7936_standard_dev_all_spring <- CO7936_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7936_standard_dev_all_spring <- CO7936_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.802644
1985 3.832749
1986 3.353474
1987 4.193416
1988 4.734483
1989 3.950378
1990 3.468374
1991 5.118588
1992 3.511918
1993 5.311126
1994 4.575242
1995 3.552697
1996 4.757735
1997 5.671268
1998 4.078121
1999 5.420418
2000 4.635885
2001 4.141790
2002 3.443972
2003 5.190409
2004 4.182243
2005 4.798876
2006 4.392900
2007 4.255123
2008 5.574087
2009 5.549686
2010 4.816305
2011 4.075459
2012 3.663627
2013 5.943932
2014 5.134123
2015 2.912585
2016 4.283511
2017 4.048996
2018 4.407727
2019 3.909838
2020 5.442965
2021 4.820355
2022 5.066184 
ggplot(CO7936_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO7936_sd_mk_spring <- mk.test(CO7936_standard_dev_all_spring$sd_2)
print(CO7936_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_spring$sd_2
## z = 0.94356, n = 39, p-value = 0.3454
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   79.0000000 6833.6666667    0.1066127
CO7936_sd_sens_spring <- sens.slope(CO7936_standard_dev_all_spring$sd_2)
print(CO7936_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_spring$sd_2
## z = 0.94356, n = 39, p-value = 0.3454
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01519521  0.03337394
## sample estimates:
## Sen's slope 
##  0.01192897

Fall temperature standard deviation

CO7936_standard_dev_all_fall <- CO7936_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7936_standard_dev_all_fall <- CO7936_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7936_standard_dev_all_fall <- CO7936_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7936_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.916629
1985 5.884137
1986 3.490781
1987 4.308378
1988 3.877165
1989 2.928018
1990 5.416582
1991 4.435321
1992 5.249867
1993 3.884032
1994 5.530286
1995 5.321321
1996 5.454961
1997 7.121094
1998 6.206828
1999 4.218235
2000 5.496021
2001 4.207251
2002 4.726965
2003 4.630696
2004 4.312028
2005 3.798870
2006 4.075571
2007 6.020582
2008 3.858402
2009 5.100863
2010 5.701597
2011 4.861906
2012 5.130162
2013 6.013378
2014 5.479466
2015 3.921206
2016 4.183132
2017 4.521044
2018 5.670414
2019 5.908079
2020 7.958004
2021 7.092743
2022 5.469556 
ggplot(CO7936_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7936 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO7936_sd_mk_fall <- mk.test(CO7936_standard_dev_all_fall$sd_2)
print(CO7936_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO7936_standard_dev_all_fall$sd_2
## z = 2.129, n = 39, p-value = 0.03325
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  177.0000000 6833.6666667    0.2388664
CO7936_sd_sens_fall <- sens.slope(CO7936_standard_dev_all_fall$sd_2)
print(CO7936_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO7936_standard_dev_all_fall$sd_2
## z = 2.129, n = 39, p-value = 0.03325
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.001631134 0.064307112
## sample estimates:
## Sen's slope 
##  0.03317649

Southern Stations

HERMIT 7 ESE CO3951

Detrending Data

#average water year temperature

CO3951_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO3951") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO3951_daily_wy_aver <- CO3951_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO3951_daily_wy_aver <- CO3951_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO3951_daily_wy_aver$aver_day_temp))

#str(CO3951_daily_wy_aver)
# try to show all years as means. 
CO3951_daily_wy_aver2 <- CO3951_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO3951_daily_wy_aver2$date_temp <- signif(CO3951_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO3951_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO3951_standard_dev <- CO3951_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO3951_standard_dev$residual)
## [1] -1.486535e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO3951_standard_dev_all <- CO3951_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3951_standard_dev_all <- CO3951_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.415412
1985 10.277513
1986 9.542602
1987 10.002142
1988 9.602689
1989 10.798841
1990 8.602389
1991 11.152040
1992 10.324649
1993 10.969684
1994 9.740607
1995 8.794524
1996 8.129266
1997 8.939283
1998 9.219493
1999 7.567821
2000 8.463930
2001 9.600737
2002 9.817103
2003 9.034250
2004 8.640154
2005 8.893249
2006 8.527566
2007 9.230474
2008 11.455350
2009 8.446360
2010 10.018255
2011 9.492496
2012 9.637339
2013 10.237010
2014 8.952825
2015 8.294359
2016 9.539250
2017 7.706123
2018 7.518636
2019 8.904952
2020 9.340397
2021 8.995937
2022 8.327062 
ggplot(CO3951_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951_sd_mk <- mk.test(CO3951_standard_dev_all$sd_2)
print(CO3951_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all$sd_2
## z = -2.5403, n = 39, p-value = 0.01107
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -211.0000000 6833.6666667   -0.2847503
CO3951_sd_sens <- sens.slope(CO3951_standard_dev_all$sd_2)
print(CO3951_sd_sens)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all$sd_2
## z = -2.5403, n = 39, p-value = 0.01107
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.063129656 -0.009152621
## sample estimates:
## Sen's slope 
## -0.03704864

Summer temperature standard deviation

CO3951_standard_dev_all_summer <- CO3951_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3951_standard_dev_all_summer <- CO3951_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.538404
1985 2.147649
1986 1.952084
1987 2.140107
1988 2.385361
1989 3.079231
1990 2.265816
1991 2.385897
1992 2.454909
1993 2.091957
1994 1.875514
1995 3.130271
1996 2.415712
1997 2.296753
1998 2.951996
1999 2.516925
2000 1.783939
2001 2.588537
2002 2.179619
2003 3.039591
2004 2.921908
2005 3.151072
2006 2.379912
2007 2.539680
2008 2.878280
2009 3.300945
2010 2.449897
2011 2.580680
2012 2.050128
2013 2.336105
2014 2.498279
2015 2.588692
2016 2.540992
2017 1.857542
2018 1.817750
2019 2.241689
2020 1.921353
2021 1.939761
2022 1.837502 
ggplot(CO3951_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO3951_sd_mk_summer <- mk.test(CO3951_standard_dev_all_summer$sd_2)
print(CO3951_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_summer$sd_2
## z = -0.67742, n = 39, p-value = 0.4981
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -57.00000000 6833.66666667   -0.07692308
CO3951_sd_sens_summer <- sens.slope(CO3951_standard_dev_all_summer$sd_2)
print(CO3951_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_summer$sd_2
## z = -0.67742, n = 39, p-value = 0.4981
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.018200664  0.007538382
## sample estimates:
##  Sen's slope 
## -0.004634582

Winter temperature standard deviation

CO3951_standard_dev_all_winter <- CO3951_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3951_standard_dev_all_winter <- CO3951_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3951_standard_dev_all_winter <- CO3951_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.967498
1985 6.536370
1986 5.105752
1987 4.663605
1988 5.470687
1989 7.132174
1990 5.469040
1991 5.619718
1992 5.388899
1993 4.067827
1994 5.674534
1995 5.172757
1996 5.241253
1997 4.668863
1998 5.126438
1999 3.945550
2000 4.904276
2001 4.578327
2002 6.002350
2003 4.653933
2004 6.098215
2005 4.823270
2006 5.185307
2007 6.222375
2008 6.918677
2009 5.365503
2010 5.785504
2011 6.284954
2012 5.283302
2013 6.712728
2014 5.121854
2015 5.531587
2016 6.076645
2017 5.967407
2018 4.496672
2019 4.774990
2020 5.511251
2021 4.821601
2022 4.714724 
ggplot(CO3951_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO3951_sd_mk_winter <- mk.test(CO3951_standard_dev_all_winter$sd_2)
print(CO3951_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_winter$sd_2
## z = -0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -27.00000000 6833.66666667   -0.03643725
CO3951_sd_sens_winter <- sens.slope(CO3951_standard_dev_all_winter$sd_2)
print(CO3951_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_winter$sd_2
## z = -0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02728328  0.02164135
## sample estimates:
##  Sen's slope 
## -0.003936932

Spring and Fall temperature standard deviation

CO3951_standard_dev_all_spring <- CO3951_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO3951_standard_dev_all_spring <- CO3951_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.551984
1985 4.046523
1986 4.475814
1987 5.155418
1988 4.514665
1989 3.266414
1990 3.040744
1991 5.869801
1992 3.868787
1993 6.069833
1994 4.496182
1995 3.971346
1996 3.583688
1997 4.393551
1998 3.897318
1999 5.021063
2000 4.592069
2001 4.829661
2002 3.045304
2003 4.746661
2004 3.283217
2005 4.774660
2006 3.611691
2007 3.818163
2008 5.608994
2009 4.793762
2010 4.809127
2011 4.323776
2012 3.690999
2013 4.085814
2014 4.566772
2015 3.015875
2016 3.954536
2017 3.829009
2018 3.077404
2019 3.415130
2020 3.502468
2021 2.731821
2022 3.692859 
ggplot(CO3951_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO3951_sd_mk_spring <- mk.test(CO3951_standard_dev_all_spring$sd_2)
print(CO3951_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_spring$sd_2
## z = -2.3952, n = 39, p-value = 0.01661
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
## -199.000000 6833.666667   -0.268556
CO3951_sd_sens_spring <- sens.slope(CO3951_standard_dev_all_spring$sd_2)
print(CO3951_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_spring$sd_2
## z = -2.3952, n = 39, p-value = 0.01661
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.059115774 -0.006797317
## sample estimates:
## Sen's slope 
## -0.03072258

Fall temperature standard deviation

CO3951_standard_dev_all_fall <- CO3951_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO3951_standard_dev_all_fall <- CO3951_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO3951_standard_dev_all_fall <- CO3951_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO3951_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.548728
1985 5.530583
1986 2.942101
1987 3.376641
1988 2.704849
1989 2.593444
1990 4.954807
1991 3.395285
1992 3.687056
1993 2.884122
1994 5.763418
1995 4.628665
1996 4.258828
1997 5.318256
1998 5.000268
1999 4.246497
2000 3.216759
2001 4.180122
2002 3.505508
2003 4.409010
2004 2.889404
2005 4.116008
2006 3.121801
2007 4.663613
2008 4.946519
2009 4.219026
2010 5.562981
2011 4.228510
2012 4.427879
2013 4.710886
2014 5.694048
2015 3.178856
2016 3.549361
2017 3.343480
2018 3.060541
2019 4.358013
2020 5.210139
2021 4.126060
2022 4.577826 
ggplot(CO3951_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO3951 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO3951_sd_mk_fall <- mk.test(CO3951_standard_dev_all_fall$sd_2)
print(CO3951_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO3951_standard_dev_all_fall$sd_2
## z = 1.0645, n = 39, p-value = 0.2871
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##   89.000000 6833.666667    0.120108
CO3951_sd_sens_fall <- sens.slope(CO3951_standard_dev_all_fall$sd_2)
print(CO3951_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO3951_standard_dev_all_fall$sd_2
## z = 1.0645, n = 39, p-value = 0.2871
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01193135  0.04435496
## sample estimates:
## Sen's slope 
##  0.01610679

PAGOSA SPRINGS CO6258

Detrending Data

#average water year temperature

CO6258_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO6258") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO6258_daily_wy_aver <- CO6258_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO6258_daily_wy_aver <- CO6258_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO6258_daily_wy_aver$aver_day_temp))

#str(CO6258_daily_wy_aver)
# try to show all years as means. 
CO6258_daily_wy_aver2 <- CO6258_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO6258_daily_wy_aver2$date_temp <- signif(CO6258_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO6258_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO6258_standard_dev <- CO6258_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO6258_standard_dev$residual)
## [1] 5.67087e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO6258_standard_dev_all <- CO6258_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO6258_standard_dev_all <- CO6258_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 9.621109
1985 8.850435
1986 8.127376
1987 8.136447
1988 9.282846
1989 9.278439
1990 9.333070
1991 9.521059
1992 8.867230
1993 8.725111
1994 9.210796
1995 7.594386
1996 8.275978
1997 8.477917
1998 9.298074
1999 7.505010
2000 8.509675
2001 9.120501
2002 9.461970
2003 8.578818
2004 8.682985
2005 8.265464
2006 8.664025
2007 9.317649
2008 9.654025
2009 8.448231
2010 9.552162
2011 9.278457
2012 9.096617
2013 10.108497
2014 8.982317
2015 8.346630
2016 9.204208
2017 8.491701
2018 8.689142
2019 9.455367
2020 6.342730 
ggplot(CO6258_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258_sd_mk <- mk.test(CO6258_standard_dev_all$sd_2)
print(CO6258_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all$sd_2
## z = 0, n = 37, p-value = 1
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##    S varS  tau 
##    0 5846    0
CO6258_sd_sens <- sens.slope(CO6258_standard_dev_all$sd_2)
print(CO6258_sd_sens)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all$sd_2
## z = 0, n = 37, p-value = 1
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01841155  0.02133537
## sample estimates:
##   Sen's slope 
## -9.500005e-05

Summer temperature standard deviation

CO6258_standard_dev_all_summer <- CO6258_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO6258_standard_dev_all_summer <- CO6258_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.892384
1985 2.715438
1986 2.591211
1987 2.509227
1988 2.800926
1989 2.890526
1990 2.720643
1991 2.759886
1992 2.859440
1993 2.634856
1994 1.609511
1995 3.033733
1996 2.114957
1997 2.291399
1998 3.466096
1999 2.339192
2000 1.556421
2001 2.215516
2002 1.798817
2003 2.795610
2004 2.271127
2005 2.848113
2006 1.810277
2007 2.411843
2008 2.442893
2009 2.719712
2010 2.249149
2011 2.347447
2012 1.751127
2013 2.083785
2014 2.157003
2015 2.092536
2016 2.370134
2017 1.874753
2018 1.971475
2019 2.536731 
ggplot(CO6258_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO6258_sd_mk_summer <- mk.test(CO6258_standard_dev_all_summer$sd_2)
print(CO6258_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_summer$sd_2
## z = -2.7923, n = 36, p-value = 0.005234
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -206.0000000 5390.0000000   -0.3269841
CO6258_sd_sens_summer <- sens.slope(CO6258_standard_dev_all_summer$sd_2)
print(CO6258_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_summer$sd_2
## z = -2.7923, n = 36, p-value = 0.005234
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.03069192 -0.00690688
## sample estimates:
## Sen's slope 
## -0.01936016

Winter temperature standard deviation

CO6258_standard_dev_all_winter <- CO6258_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO6258_standard_dev_all_winter <- CO6258_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO6258_standard_dev_all_winter <- CO6258_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.539382
1985 5.221322
1986 5.151056
1987 4.006869
1988 5.463126
1989 6.081710
1990 5.501666
1991 5.959856
1992 4.854584
1993 4.839559
1994 4.459638
1995 4.086126
1996 4.602684
1997 4.824240
1998 4.287481
1999 3.999260
2000 4.513163
2001 3.948779
2002 5.199542
2003 3.698090
2004 5.469668
2005 3.678975
2006 4.423239
2007 5.636067
2008 6.071441
2009 4.686136
2010 4.810118
2011 5.553928
2012 4.251732
2013 6.212543
2014 4.618260
2015 4.988689
2016 5.296415
2017 5.513993
2018 4.031030
2019 4.772180
2020 4.874283 
ggplot(CO6258_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO6258_sd_mk_winter <- mk.test(CO6258_standard_dev_all_winter$sd_2)
print(CO6258_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_winter$sd_2
## z = -0.48392, n = 37, p-value = 0.6284
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -38.00000000 5846.00000000   -0.05705706
CO6258_sd_sens_winter <- sens.slope(CO6258_standard_dev_all_winter$sd_2)
print(CO6258_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_winter$sd_2
## z = -0.48392, n = 37, p-value = 0.6284
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.03081031  0.01926972
## sample estimates:
##  Sen's slope 
## -0.005776429

Spring and Fall temperature standard deviation

CO6258_standard_dev_all_spring <- CO6258_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO6258_standard_dev_all_spring <- CO6258_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.490422
1985 3.431979
1986 3.874644
1987 3.460695
1988 4.169099
1989 3.249797
1990 2.531130
1991 4.202630
1992 2.908047
1993 4.142813
1994 4.202854
1995 3.393549
1996 4.035316
1997 4.329517
1998 4.058368
1999 4.509461
2000 4.104678
2001 4.193235
2002 3.070249
2003 4.649083
2004 3.374868
2005 4.328198
2006 3.599633
2007 3.695058
2008 4.056428
2009 4.644815
2010 4.071791
2011 3.648932
2012 3.716872
2013 4.597192
2014 4.657398
2015 3.154261
2016 3.833665
2017 3.935358
2018 4.322978
2019 3.446913 
ggplot(CO6258_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO6258_sd_mk_spring <- mk.test(CO6258_standard_dev_all_spring$sd_2)
print(CO6258_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_spring$sd_2
## z = 0.64018, n = 36, p-value = 0.5221
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.800000e+01 5.390000e+03 7.619048e-02
CO6258_sd_sens_spring <- sens.slope(CO6258_standard_dev_all_spring$sd_2)
print(CO6258_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_spring$sd_2
## z = 0.64018, n = 36, p-value = 0.5221
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01196439  0.02627478
## sample estimates:
## Sen's slope 
## 0.007762508

Fall temperature standard deviation

CO6258_standard_dev_all_fall <- CO6258_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO6258_standard_dev_all_fall <- CO6258_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO6258_standard_dev_all_fall <- CO6258_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO6258_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.606103
1985 5.482743
1986 3.369019
1987 3.990539
1988 3.365853
1989 2.964852
1990 5.488376
1991 3.860703
1992 3.882496
1993 3.415086
1994 4.706409
1995 4.284017
1996 4.226385
1997 5.675789
1998 6.103044
1999 3.435030
2000 4.168726
2001 4.217347
2002 3.702107
2003 3.997812
2004 3.321224
2005 4.452367
2006 3.323123
2007 4.953457
2008 3.553112
2009 4.135746
2010 5.534732
2011 4.032637
2012 4.617749
2013 5.070778
2014 5.318232
2015 3.571338
2016 3.472990
2017 3.891573
2018 4.668249
2019 4.970637
2020 5.544065 
ggplot(CO6258_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO6258 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO6258_sd_mk_fall <- mk.test(CO6258_standard_dev_all_fall$sd_2)
print(CO6258_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO6258_standard_dev_all_fall$sd_2
## z = 1.5825, n = 37, p-value = 0.1135
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  122.0000000 5846.0000000    0.1831832
CO6258_sd_sens_fall <- sens.slope(CO6258_standard_dev_all_fall$sd_2)
print(CO6258_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO6258_standard_dev_all_fall$sd_2
## z = 1.5825, n = 37, p-value = 0.1135
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.005727141  0.047076532
## sample estimates:
## Sen's slope 
##  0.02035048

SILVERTON CO7656

Detrending Data

#average water year temperature

CO7656_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO7656") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

CO7656_temp_xts <- xts(CO7656_yearly_wy_aver$avg_T_c, order.by = CO7656_yearly_wy_aver$Date)

dygraph(CO7656_temp_xts) %>%
  dyAxis("y", label = "Daily temperature (°C)") 
#Average temperature by day for all water years:

CO7656_daily_wy_aver <- CO7656_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO7656_daily_wy_aver <- CO7656_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO7656_daily_wy_aver$aver_day_temp))

#str(CO7656_daily_wy_aver)
# try to show all years as means. 
CO7656_daily_wy_aver2 <- CO7656_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO7656_daily_wy_aver2$date_temp <- signif(CO7656_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO7656_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO7656_standard_dev <- CO7656_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO7656_standard_dev$residual)
## [1] -3.691746e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO7656_standard_dev_all <- CO7656_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7656_standard_dev_all <- CO7656_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 9.800356
1985 9.286243
1986 8.190448
1987 8.552234
1988 9.163062
1989 9.244891
1990 8.705592
1991 9.568655
1992 8.473835
1993 8.837258
1994 9.549772
1995 8.757576
1996 8.542260
1997 9.124829
1998 9.393025
1999 7.785032
2000 8.493772
2001 9.124296
2002 9.260411
2003 8.916683
2004 8.666508
2005 8.031985
2006 8.595742
2007 9.260942
2008 9.602807
2009 8.642532
2010 9.778187
2011 9.699779
2012 10.200554
2013 10.042336
2014 9.013978
2015 8.191825
2016 9.335647
2017 8.888276
2018 9.205000
2019 11.507607
2020 9.458099
2021 9.407391
2022 9.248475 
ggplot(CO7656_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656_sd_mk <- mk.test(CO7656_standard_dev_all$sd_2)
print(CO7656_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all$sd_2
## z = 1.6452, n = 39, p-value = 0.09993
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  137.0000000 6833.6666667    0.1848853
CO7656_sd_sens <- sens.slope(CO7656_standard_dev_all$sd_2)
print(CO7656_sd_sens)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all$sd_2
## z = 1.6452, n = 39, p-value = 0.09993
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.00318769  0.03233149
## sample estimates:
## Sen's slope 
##  0.01553101

Summer temperature standard deviation

CO7656_standard_dev_all_summer <- CO7656_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7656_standard_dev_all_summer <- CO7656_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.868027
1985 2.392375
1986 2.106835
1987 2.443995
1988 3.004058
1989 2.703892
1990 2.400346
1991 2.388603
1992 2.640532
1993 2.633475
1994 1.963849
1995 3.130228
1996 2.187039
1997 2.342750
1998 3.061265
1999 2.716929
2000 1.875917
2001 2.472874
2002 2.269983
2003 2.452708
2004 2.134423
2005 2.633342
2006 3.584119
2007 2.688492
2008 2.552210
2009 2.690311
2010 2.333727
2011 2.388600
2012 2.682441
2013 2.110861
2014 2.419230
2015 2.100976
2016 2.183482
2017 1.958013
2018 4.401762
2019 5.124136
2020 2.310070
2021 2.053687
2022 2.089330 
ggplot(CO7656_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO7656_sd_mk_summer <- mk.test(CO7656_standard_dev_all_summer$sd_2)
print(CO7656_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_summer$sd_2
## z = -1.3548, n = 39, p-value = 0.1755
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## -113.0000000 6833.6666667   -0.1524966
CO7656_sd_sens_summer <- sens.slope(CO7656_standard_dev_all_summer$sd_2)
print(CO7656_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_summer$sd_2
## z = -1.3548, n = 39, p-value = 0.1755
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.017626864  0.004829158
## sample estimates:
##  Sen's slope 
## -0.007146811

Winter temperature standard deviation

CO7656_standard_dev_all_winter <- CO7656_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7656_standard_dev_all_winter <- CO7656_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7656_standard_dev_all_winter <- CO7656_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.825255
1985 5.791930
1986 5.050040
1987 4.148185
1988 5.340055
1989 6.451878
1990 5.408032
1991 5.549401
1992 4.346191
1993 4.615703
1994 4.994627
1995 5.023145
1996 5.152052
1997 5.188795
1998 5.001190
1999 4.579430
2000 5.007575
2001 4.586201
2002 5.744089
2003 4.401514
2004 5.709564
2005 4.848385
2006 4.140175
2007 5.949719
2008 6.572679
2009 5.229943
2010 5.141441
2011 6.174123
2012 4.602861
2013 6.207070
2014 4.883980
2015 5.137885
2016 5.442594
2017 6.400228
2018 5.314075
2019 5.040701
2020 5.359974
2021 4.899933
2022 5.611303 
ggplot(CO7656_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO7656_sd_mk_winter <- mk.test(CO7656_standard_dev_all_winter$sd_2)
print(CO7656_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_winter$sd_2
## z = 0.58065, n = 39, p-value = 0.5615
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 4.900000e+01 6.833667e+03 6.612686e-02
CO7656_sd_sens_winter <- sens.slope(CO7656_standard_dev_all_winter$sd_2)
print(CO7656_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_winter$sd_2
## z = 0.58065, n = 39, p-value = 0.5615
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01262523  0.02718868
## sample estimates:
## Sen's slope 
## 0.005991615

Spring and Fall temperature standard deviation

CO7656_standard_dev_all_spring <- CO7656_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO7656_standard_dev_all_spring <- CO7656_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.054476
1985 4.008601
1986 4.073680
1987 4.321120
1988 4.256750
1989 4.060835
1990 2.824572
1991 5.168149
1992 3.240501
1993 4.642508
1994 4.425972
1995 3.964639
1996 5.097443
1997 5.343814
1998 4.994585
1999 5.109918
2000 4.541188
2001 5.085421
2002 2.797888
2003 5.391350
2004 3.317837
2005 4.333056
2006 3.763982
2007 4.179142
2008 4.768617
2009 5.606861
2010 4.816244
2011 4.474883
2012 5.612087
2013 4.344215
2014 5.190846
2015 3.112420
2016 3.943439
2017 4.039696
2018 3.633129
2019 5.409801
2020 4.389195
2021 3.446429
2022 4.412791 
ggplot(CO7656_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO7656_sd_mk_spring <- mk.test(CO7656_standard_dev_all_spring$sd_2)
print(CO7656_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_spring$sd_2
## z = 0, n = 39, p-value = 1
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 1.000000e+00 6.833667e+03 1.349528e-03
CO7656_sd_sens_spring <- sens.slope(CO7656_standard_dev_all_spring$sd_2)
print(CO7656_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_spring$sd_2
## z = 0, n = 39, p-value = 1
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02669756  0.02437292
## sample estimates:
##  Sen's slope 
## 0.0006631563

Fall temperature standard deviation

CO7656_standard_dev_all_fall <- CO7656_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO7656_standard_dev_all_fall <- CO7656_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO7656_standard_dev_all_fall <- CO7656_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO7656_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.022752
1985 5.478970
1986 3.383785
1987 3.633847
1988 3.061820
1989 2.803181
1990 5.088947
1991 4.359722
1992 4.102061
1993 2.792590
1994 4.468378
1995 4.661300
1996 4.053867
1997 6.027621
1998 5.949864
1999 3.629842
2000 3.648252
2001 4.173251
2002 3.529245
2003 4.400172
2004 3.212547
2005 3.843181
2006 3.541248
2007 4.960946
2008 3.902991
2009 3.800532
2010 5.489741
2011 4.209702
2012 4.296304
2013 4.711007
2014 5.593147
2015 2.989476
2016 3.268318
2017 3.529796
2018 3.449304
2019 4.584459
2020 5.780197
2021 3.692796
2022 4.965113 
ggplot(CO7656_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO7656 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO7656_sd_mk_fall <- mk.test(CO7656_standard_dev_all_fall$sd_2)
print(CO7656_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO7656_standard_dev_all_fall$sd_2
## z = 1.0887, n = 39, p-value = 0.2763
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
##   91.000000 6833.666667    0.122807
CO7656_sd_sens_fall <- sens.slope(CO7656_standard_dev_all_fall$sd_2)
print(CO7656_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO7656_standard_dev_all_fall$sd_2
## z = 1.0887, n = 39, p-value = 0.2763
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01411259  0.04886932
## sample estimates:
## Sen's slope 
##  0.01722729

VALLECITO DAM CO8582

Detrending Data

#average water year temperature

CO8582_yearly_wy_aver <- COOP_stations_clean %>% 
  filter(station == "CO8582") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))

CO8582_temp_xts <- xts(CO8582_yearly_wy_aver$avg_T_c, order.by = CO8582_yearly_wy_aver$Date)

dygraph(CO8582_temp_xts) %>%
  dyAxis("y", label = "Daily temperature (°C)") 
#Average temperature by day for all water years:

CO8582_daily_wy_aver <- CO8582_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO8582_daily_wy_aver <- CO8582_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO8582_daily_wy_aver$aver_day_temp))

#str(CO8582_daily_wy_aver)
# try to show all years as means. 
CO8582_daily_wy_aver2 <- CO8582_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO8582_daily_wy_aver2$date_temp <- signif(CO8582_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO8582_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO8582_standard_dev <- CO8582_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO8582_standard_dev$residual)
## [1] -6.23937e-16

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO8582_standard_dev_all <- CO8582_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8582_standard_dev_all <- CO8582_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 9.383586
1985 8.747479
1986 7.843649
1987 8.290687
1988 9.157787
1989 9.113597
1990 8.535844
1991 9.094549
1992 8.416527
1993 8.647521
1994 9.052431
1995 8.578530
1996 8.579043
1997 9.107490
1998 9.383756
1999 7.505741
2000 8.676708
2001 9.337435
2002 9.637440
2003 8.889594
2004 9.112936
2005 8.538574
2006 8.461740
2007 9.192161
2008 9.868298
2009 8.660182
2010 9.849250
2011 9.352170
2012 9.066766
2013 9.972813
2014 8.277648
2015 7.734352
2016 8.860104
2017 8.171528
2018 8.242078
2019 9.613219
2020 9.334744
2021 9.524267
2022 8.975630 
ggplot(CO8582_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582_sd_mk <- mk.test(CO8582_standard_dev_all$sd_2)
print(CO8582_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all$sd_2
## z = 0.94356, n = 39, p-value = 0.3454
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##   79.0000000 6833.6666667    0.1066127
CO8582_sd_sens <- sens.slope(CO8582_standard_dev_all$sd_2)
print(CO8582_sd_sens)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all$sd_2
## z = 0.94356, n = 39, p-value = 0.3454
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.009377773  0.030402032
## sample estimates:
## Sen's slope 
##  0.01084421

Summer temperature standard deviation

CO8582_standard_dev_all_summer <- CO8582_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8582_standard_dev_all_summer <- CO8582_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.813935
1985 2.637505
1986 2.404488
1987 2.448139
1988 2.532795
1989 2.833055
1990 2.690463
1991 2.554909
1992 2.738031
1993 2.938426
1994 1.921098
1995 3.603464
1996 1.873828
1997 2.502794
1998 3.417954
1999 2.629655
2000 1.999305
2001 2.422970
2002 2.118775
2003 2.967155
2004 2.537188
2005 3.173025
2006 2.199192
2007 2.735728
2008 2.789759
2009 3.198484
2010 2.485891
2011 2.499998
2012 1.885593
2013 2.159771
2014 2.354197
2015 2.552712
2016 2.537342
2017 2.201989
2018 2.089558
2019 2.877991
2020 2.368652
2021 2.142237
2022 2.401416 
ggplot(CO8582_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO8582_sd_mk_summer <- mk.test(CO8582_standard_dev_all_summer$sd_2)
print(CO8582_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_summer$sd_2
## z = -1.5, n = 39, p-value = 0.1336
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##           S        varS         tau 
## -125.000000 6833.666667   -0.168691
CO8582_sd_sens_summer <- sens.slope(CO8582_standard_dev_all_summer$sd_2)
print(CO8582_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_summer$sd_2
## z = -1.5, n = 39, p-value = 0.1336
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.018363877  0.003232946
## sample estimates:
##  Sen's slope 
## -0.008549483

Winter temperature standard deviation

CO8582_standard_dev_all_winter <- CO8582_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8582_standard_dev_all_winter <- CO8582_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8582_standard_dev_all_winter <- CO8582_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.310241
1985 5.045484
1986 4.862671
1987 4.267177
1988 5.301464
1989 6.017133
1990 5.570090
1991 5.607337
1992 4.574863
1993 4.345881
1994 4.847218
1995 4.719710
1996 5.023067
1997 5.487041
1998 4.711243
1999 3.973766
2000 4.915390
2001 4.399419
2002 5.274510
2003 4.056802
2004 6.212342
2005 4.020457
2006 4.485868
2007 5.860868
2008 6.447097
2009 5.180071
2010 5.170040
2011 5.570027
2012 4.282039
2013 6.514920
2014 4.645525
2015 4.899033
2016 5.148721
2017 5.486057
2018 4.427436
2019 4.891143
2020 4.681322
2021 5.170427
2022 5.084264 
ggplot(CO8582_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO8582_sd_mk_winter <- mk.test(CO8582_standard_dev_all_winter$sd_2)
print(CO8582_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_winter$sd_2
## z = -0.048387, n = 39, p-value = 0.9614
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
## -5.000000e+00  6.833667e+03 -6.747638e-03
CO8582_sd_sens_winter <- sens.slope(CO8582_standard_dev_all_winter$sd_2)
print(CO8582_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_winter$sd_2
## z = -0.048387, n = 39, p-value = 0.9614
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01834464  0.02017638
## sample estimates:
##   Sen's slope 
## -0.0008037351

Spring and Fall temperature standard deviation

CO8582_standard_dev_all_spring <- CO8582_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO8582_standard_dev_all_spring <- CO8582_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 5.582789
1985 3.015812
1986 4.116645
1987 3.458456
1988 4.339430
1989 3.524440
1990 2.697893
1991 4.018877
1992 3.190276
1993 4.111149
1994 4.447072
1995 3.310654
1996 4.515146
1997 5.105015
1998 4.653657
1999 5.192960
2000 4.346563
2001 4.998349
2002 3.362700
2003 4.994796
2004 3.667480
2005 4.509563
2006 3.891584
2007 4.069403
2008 4.217933
2009 5.035479
2010 4.441639
2011 3.826997
2012 3.535582
2013 4.457302
2014 4.370804
2015 2.978738
2016 3.639620
2017 3.420019
2018 3.950565
2019 3.654574
2020 3.915741
2021 3.462782
2022 4.020531 
ggplot(CO8582_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO8582_sd_mk_spring <- mk.test(CO8582_standard_dev_all_spring$sd_2)
print(CO8582_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_spring$sd_2
## z = -0.50807, n = 39, p-value = 0.6114
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##             S          varS           tau 
##  -43.00000000 6833.66666667   -0.05802969
CO8582_sd_sens_spring <- sens.slope(CO8582_standard_dev_all_spring$sd_2)
print(CO8582_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_spring$sd_2
## z = -0.50807, n = 39, p-value = 0.6114
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02656507  0.01642275
## sample estimates:
##  Sen's slope 
## -0.005794715

Fall temperature standard deviation

CO8582_standard_dev_all_fall <- CO8582_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO8582_standard_dev_all_fall <- CO8582_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO8582_standard_dev_all_fall <- CO8582_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO8582_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.669659
1985 5.166616
1986 3.385639
1987 4.093466
1988 3.456959
1989 3.196028
1990 4.089225
1991 3.758301
1992 4.580766
1993 3.126566
1994 4.895095
1995 5.891050
1996 4.100669
1997 6.316312
1998 5.840303
1999 3.432794
2000 4.329076
2001 4.501070
2002 4.167554
2003 4.355277
2004 3.807048
2005 4.477278
2006 3.337021
2007 5.088793
2008 3.811342
2009 4.364315
2010 5.765552
2011 4.211652
2012 4.734032
2013 4.684849
2014 5.326233
2015 3.908040
2016 3.040585
2017 4.058289
2018 3.933025
2019 5.210893
2020 6.345771
2021 4.797724
2022 5.607850 
ggplot(CO8582_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO8582 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO8582_sd_mk_fall <- mk.test(CO8582_standard_dev_all_fall$sd_2)
print(CO8582_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO8582_standard_dev_all_fall$sd_2
## z = 1.9113, n = 39, p-value = 0.05597
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  159.0000000 6833.6666667    0.2145749
CO8582_sd_sens_fall <- sens.slope(CO8582_standard_dev_all_fall$sd_2)
print(CO8582_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO8582_standard_dev_all_fall$sd_2
## z = 1.9113, n = 39, p-value = 0.05597
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.001413698  0.049766343
## sample estimates:
## Sen's slope 
##  0.02458039

DEL NORTE 2E CO2184

Detrending Data

#average water year temperature

CO2184_yearly_wy_aver <- COOP_stations_clean %>%
  filter(station == "CO2184") %>% 
  group_by(waterYear) %>%
  filter(waterYear >= 1984) %>% 
  mutate(aver_ann_temp = mean(avg_T_c))
#Average temperature by day for all water years:

CO2184_daily_wy_aver <- CO2184_yearly_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(aver_day_temp = mean(aver_ann_temp))

#average mean temperature by day for the period of record:

CO2184_daily_wy_aver <- CO2184_daily_wy_aver %>% 
  group_by(daymonth) %>% 
  mutate(all_ave_temp = mean(CO2184_daily_wy_aver$aver_day_temp))

#str(CO2184_daily_wy_aver)
# try to show all years as means. 
CO2184_daily_wy_aver2 <- CO2184_daily_wy_aver %>% 
  #filter(waterYear == "1987" | waterYear == "2021") %>%
  group_by(waterDay) %>%
  mutate(date_temp = mean(avg_T_c))
  
CO2184_daily_wy_aver2$date_temp <- signif(CO2184_daily_wy_aver2$date_temp,3) #reduce the sig figs

ggplot(CO2184_daily_wy_aver2, aes(x = waterDay, y = date_temp))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  #geom_smooth(method = "lm", se=FALSE) +
  ylab('Average Daily temperature (°C)') + 
  xlab('Day of water year')

Standard deviation- Determining residuals

CO2184_standard_dev <- CO2184_daily_wy_aver %>% 
  group_by(waterYear) %>% 
  mutate(residual = (all_ave_temp-aver_ann_temp)+avg_T_c-aver_day_temp) %>% 
  mutate(deviation = abs(residual-lag(residual)))

mean(CO2184_standard_dev$residual)
## [1] -3.33898e-17

The mean of the residuals is close enough to zero

Calculating standard deviation for the timeseries

CO2184_standard_dev_all <- CO2184_standard_dev %>% 
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO2184_standard_dev_all <- CO2184_standard_dev_all %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 10.851354
1985 9.322284
1986 8.117289
1987 8.930130
1988 10.122162
1989 9.678175
1990 8.420811
1991 9.283741
1992 10.283845
1993 9.683157
1994 9.228863
1995 7.734984
1996 8.180139
1997 8.928496
1998 9.265695
1999 7.731617
2000 8.352092
2001 9.255820
2002 9.754404
2003 8.884745
2004 9.145182
2005 8.664351
2006 8.681983
2007 10.102809
2008 10.762096
2009 8.982543
2010 9.892880
2011 9.413105
2012 10.074149
2013 11.257504
2014 10.017296
2015 8.434601
2016 8.996848
2017 8.840504
2018 8.044550
2019 9.547470
2020 10.455481
2021 9.952095
2022 9.426592 
#CALLING THIS something different
CO2184_all_V2 <- ggplot(CO2184_standard_dev_all, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')


print(CO2184_all_V2)

CO2184_sd_mk <- mk.test(CO2184_standard_dev_all$sd_2)
print(CO2184_sd_mk)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all$sd_2
## z = 0.72581, n = 39, p-value = 0.468
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 6.100000e+01 6.833667e+03 8.232119e-02
CO2184_sd_sens <- sens.slope(CO2184_standard_dev_all$sd_2)
print(CO2184_sd_sens)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all$sd_2
## z = 0.72581, n = 39, p-value = 0.468
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01500289  0.03885031
## sample estimates:
## Sen's slope 
##  0.01041623

Summer temperature standard deviation

CO2184_standard_dev_all_summer <- CO2184_standard_dev %>%
  filter(waterDay >= 244 & waterDay <= 335) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO2184_standard_dev_all_summer <- CO2184_standard_dev_all_summer %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_summer %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 2.391702
1985 2.143789
1986 1.884837
1987 1.827067
1988 2.187482
1989 2.292524
1990 2.135460
1991 2.107556
1992 2.327857
1993 2.347764
1994 1.627795
1995 2.551853
1996 2.100201
1997 1.926402
1998 2.600348
1999 2.410172
2000 1.600572
2001 2.019630
2002 1.845975
2003 2.801925
2004 2.332780
2005 2.619140
2006 1.996965
2007 2.363320
2008 2.356730
2009 2.488327
2010 1.954418
2011 2.097681
2012 1.719314
2013 2.269814
2014 2.174603
2015 2.003585
2016 2.357237
2017 2.187614
2018 2.098239
2019 2.398375
2020 2.058420
2021 2.540783
2022 2.216535 
ggplot(CO2184_standard_dev_all_summer, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Summer standard deviations.

CO2184_sd_mk_summer <- mk.test(CO2184_standard_dev_all_summer$sd_2)
print(CO2184_sd_mk_summer)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_summer$sd_2
## z = 0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 5.300000e+01 6.833667e+03 7.152497e-02
CO2184_sd_sens_summer <- sens.slope(CO2184_standard_dev_all_summer$sd_2)
print(CO2184_sd_sens_summer)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_summer$sd_2
## z = 0.62904, n = 39, p-value = 0.5293
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.005919914  0.010498988
## sample estimates:
## Sen's slope 
## 0.002611781

Winter temperature standard deviation

CO2184_standard_dev_all_winter <- CO2184_standard_dev %>%
  filter(waterDay >= 32 & waterDay <= 182) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO2184_standard_dev_all_winter <- CO2184_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO2184_standard_dev_all_winter <- CO2184_standard_dev_all_winter %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_winter %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.386911
1985 5.608819
1986 5.599562
1987 4.291683
1988 6.239150
1989 6.819955
1990 5.352755
1991 6.094841
1992 5.262553
1993 5.601879
1994 5.076673
1995 4.452943
1996 4.844600
1997 5.565064
1998 4.839934
1999 4.355520
2000 4.682674
2001 4.573170
2002 5.564876
2003 4.760996
2004 6.416285
2005 4.856311
2006 4.991127
2007 7.036398
2008 7.257602
2009 5.728112
2010 5.137456
2011 5.784141
2012 6.033598
2013 7.479222
2014 6.422942
2015 5.408038
2016 5.599011
2017 6.234240
2018 3.817411
2019 5.154514
2020 6.145269
2021 5.186261
2022 5.479436 
ggplot(CO2184_standard_dev_all_winter, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Winter standard deviations.

CO2184_sd_mk_winter <- mk.test(CO2184_standard_dev_all_winter$sd_2)
print(CO2184_sd_mk_winter)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_winter$sd_2
## z = 0.41129, n = 39, p-value = 0.6809
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 3.500000e+01 6.833667e+03 4.723347e-02
CO2184_sd_sens_winter <- sens.slope(CO2184_standard_dev_all_winter$sd_2)
print(CO2184_sd_sens_winter)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_winter$sd_2
## z = 0.41129, n = 39, p-value = 0.6809
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.02091837  0.03094743
## sample estimates:
## Sen's slope 
## 0.005319822

Spring and Fall temperature standard deviation

CO2184_standard_dev_all_spring <- CO2184_standard_dev %>%
  filter(waterDay >= 183 & waterDay <= 243) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

CO2184_standard_dev_all_spring <- CO2184_standard_dev_all_spring %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_spring %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 6.010688
1985 3.158078
1986 3.580804
1987 3.721681
1988 3.690954
1989 4.098044
1990 3.175753
1991 3.959011
1992 3.160167
1993 3.877660
1994 4.477492
1995 3.650002
1996 5.013044
1997 5.023724
1998 4.271152
1999 4.880544
2000 4.986695
2001 4.063152
2002 2.903261
2003 4.981343
2004 4.166573
2005 4.556008
2006 4.026164
2007 4.280817
2008 4.410486
2009 4.732733
2010 4.434111
2011 4.067431
2012 3.576749
2013 4.988453
2014 4.558860
2015 2.766850
2016 3.960783
2017 3.949906
2018 3.460995
2019 3.452481
2020 4.130015
2021 3.997066
2022 4.589323 
ggplot(CO2184_standard_dev_all_spring, aes(x = waterYear, y = sd_2))+
  geom_line(size= 0.7) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Jun-Aug standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Spring standard deviations.

CO2184_sd_mk_spring <- mk.test(CO2184_standard_dev_all_spring$sd_2)
print(CO2184_sd_mk_spring)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_spring$sd_2
## z = 0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
## 2.700000e+01 6.833667e+03 3.643725e-02
CO2184_sd_sens_spring <- sens.slope(CO2184_standard_dev_all_spring$sd_2)
print(CO2184_sd_sens_spring)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_spring$sd_2
## z = 0.31452, n = 39, p-value = 0.7531
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  -0.01775251  0.02589369
## sample estimates:
## Sen's slope 
##  0.00361402

Fall temperature standard deviation

CO2184_standard_dev_all_fall <- CO2184_standard_dev %>%
  filter(waterDay >= 336 | waterDay <= 31) %>% # this might be better off as daymonth rather than day of water year due to leap year
  group_by(waterYear) %>% 
  mutate(nmbr = n())

# Nope. This did some weird stuff with twice the observations.
# CO2184_standard_dev_all_fall <- CO2184_standard_dev %>%
#   filter(daymonth >= "01-11" & daymonth <= "31-03") %>%
#   group_by(waterYear) %>% 
#   mutate(nmbr = n())


CO2184_standard_dev_all_fall <- CO2184_standard_dev_all_fall %>% 
  group_by(waterYear) %>% 
  mutate(resid_mean = mean(residual)) %>%
  mutate(sd_1 = residual-resid_mean) %>% 
  mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
  distinct(sd_2, .keep_all = TRUE) %>% 
   select(waterYear, sd_2)

CO2184_standard_dev_all_fall %>% 
  kable(.,'html') %>%
  kable_styling() %>%
  scroll_box(width='250px',height='500px')
waterYear sd_2
1984 3.514378
1985 5.100490
1986 2.876471
1987 3.541120
1988 2.911261
1989 2.511399
1990 5.108661
1991 3.739867
1992 4.316150
1993 2.494880
1994 5.563351
1995 4.320884
1996 3.875344
1997 5.921244
1998 5.100024
1999 3.938286
2000 4.327829
2001 4.416964
2002 4.134282
2003 3.962053
2004 3.370771
2005 4.537013
2006 3.407062
2007 5.065168
2008 3.990375
2009 4.177926
2010 5.993435
2011 4.696400
2012 4.590017
2013 5.248431
2014 5.066031
2015 3.685121
2016 3.383254
2017 3.779273
2018 4.585144
2019 4.670623
2020 6.095844
2021 5.906143
2022 5.443732 
ggplot(CO2184_standard_dev_all_fall, aes(x = waterYear, y = sd_2))+#, color = waterYear)) +
  geom_line(size= 0.7) +
  #geom_line(aes) +
  theme_few() +
  geom_smooth(method = "lm", se=FALSE) +
  ylab('SD') + 
  xlab('Water year')

CO2184 Nov-Mar standard deviation for water years 1984-2022

Mann-Kendall & Sen’s Slope

Fall standard deviations.

CO2184_sd_mk_fall <- mk.test(CO2184_standard_dev_all_fall$sd_2)
print(CO2184_sd_mk_fall)
## 
##  Mann-Kendall trend test
## 
## data:  CO2184_standard_dev_all_fall$sd_2
## z = 2.5645, n = 39, p-value = 0.01033
## alternative hypothesis: true S is not equal to 0
## sample estimates:
##            S         varS          tau 
##  213.0000000 6833.6666667    0.2874494
CO2184_sd_sens_fall <- sens.slope(CO2184_standard_dev_all_fall$sd_2)
print(CO2184_sd_sens_fall)
## 
##  Sen's slope
## 
## data:  CO2184_standard_dev_all_fall$sd_2
## z = 2.5645, n = 39, p-value = 0.01033
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
##  0.009585968 0.062560256
## sample estimates:
## Sen's slope 
##  0.03404516